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Abstract 

Several applications of spectral methods to problems re- 
lated to the relativistic astrophysics of compact objects 
are presented. Based on a proper definition of the an- 
alytical properties of regular tensorial functions we have 
developed a spectral method in a general sphericallike co- 
ordinate system. The applications include the investiga- 
tion of spherically symmetric neutron star collapse as well 
as the solution of the coupled 2D-Einstein-Maxwell equa- 
tions for magnetized, rapidly rotating neutron stars. In 
both cases the resulting codes are efficient and give results 
typically several orders of magnitude more accurate than 
equivalent codes based on finite difference schemes. We 
further report the current status of a 3D-code aiming at 
the simulation of non-axisymmetric neutron star collapse 
where we have chosen a tensor based numerical scheme. 

Key words: spectral methods, numerical relativity, neu- 
tron stars, black holes, gravitational collapse. 
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1 Introduction 

Compact objects in astrophysics such as neutron stars 
and black holes are subjected to the strong field regime 
of gravitation and have hence to be treated within the 
framework of general relativity. The growing interest in 
the numerical solution of the Einstein equations for astro- 
physically relevant systems has given rise to a new branch 
of computational physics — numerical relativity [[l], |^, ^| . 
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This development is due to the increasingly powerful com- 
putational resources which make these problems accessi- 
ble to a numerical investigation. It is further stimulated 
by the prospects of gravitational wave astronomy which 
will turn into an observational science toward the end of 
this decade thanks to gravitational wave observatories like 
LIGO, VIRGO and GEO600 that are now under construe- 
tion g||. 

We use the (3+1) -formalism of general relativity |jj which 
consists in foliating spacetime into a sequence of space- 
like hypersurfaces which represent curved three-space at 
a fixed coordinate time t. The fabric of spacetime is then 
determined by the three-metric hij and four additional 
quantities, the lapse function N and the shift vector N l 
which fix the propagation of the spacelike hypersurfaces 
in time and the change of the spatial coordinate system 
between adjacent hypersurfaces. This Hamilton type ap- 
proach to general relativity results in a temporal first order 
evolution scheme for the dynamical variables which is com- 
pleted by some constraint equations which ensure the con- 
sistency of gravitational and matter fields. Furthermore N 
and N l have to be determined by the choice of appropiate 
gauge conditions which typically lead to elliptic equations 
that have to be solved at each time step. For stationary 
configurations all time derivatives vanish and one obtains 
a system of coupled elliptic equations for the gravitational 
fields. The efficient solution of elliptic equations is hence 
of central interest for us. 

Let us consider a covariant Poisson equation N' i u = S in a 
conformally flat axisymmetric space where the line element 
reads 



(1) 



dl 2 =A 4 (r, 9) (dr 2 +r 2 d9 2 +r 2 sin 2 



The former equation can be rewritten to yield a Poisson- 
like equation for N where we have isolated the flat space 
Laplacian Af and contributed the curvature terms to the 
source. Here a denotes In A. 

(2) A f N = S with S = A i S-2(d r ad r N+\dead e N). 
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This equation has to be solved by iteration. The solution 
of Af N — S at each iteration has hence to be accomplished 
sufficiently fast in order to keep the total computation cost 
at a reasonable level. 

After outlining the basic features of our spectral method 
[||, [|, we will proceed in a first step to the investigation 
of black hole formation due to spherically symmetric neu- 
tron star collapse which has proved the high aptitude of 
spectral methods in this field Jl(| |ll], [l2|. The second 
part is devoted to the study of axisymmetric stationary 
rotating bodies which has been applied to model rapidly 
rotating neutron stars [^3[ 0. This work has been ex- 
tended recently to include strong magnetic fields for the 
first time into neutron star models [l5| |. Special emphasis 
in all cases has been put on the extensive use of external 
and intrinsic tests 



which can be written in terms of (r, 9, <f>) as 



liet 17 



|18[ of the self-consistency and 
the attained accuracy of the numerical results. The result- 
ing neutron star models provide us with the required initial 
value models for the investigation of 3D-gravitational col- 
lapse of neutron stars which will reveal the whole range 
of gravitational wave emission associated with this phe- 
nomenon. We give an overview about the inset of spectral 
methods in this project which is currently in work. Here 
a new method for the efficient inversion of a generalized 
3D-vector Poisson equation is a first major result. 

2 Spectral methods in general rel- 
ativity 

2.1 Coordinates and regularity conditions 

The spacelike hypersurfaces stemming from the former 
choice of the (3+1) -formalism are conceived to describe 
some asymptotically flat space, containing a compact, 
mostly starlike object. The natural choice is thus a sphe- 
ricallike coordinate system (r, 0, 4>). 

The pseudosingularities which appear in this case can be 
overcome by a proper definition of regularity conditions 
of tensorial quantities. A consequent application of parity 
rules derived from these conditions allows further to opti- 
mize code efficiency and precision. 

We consider the related Cartesian type coordinate system 
(x,y,z) = (r sin 9 cos </>, r sin 9 sin 4>, r cos 9) . We define a 
tensorial quantity T^...^ in spherical coordinates (r,9,<fi) 
to be regular, if its components fi l ...i N with respect to 
Cartesian coordinates (x, y, z) are regular in the sense that 
they can be expanded into a polynomial sum of the type 



N 



(3) 



fii...i N {x, y, z) 



i,j : k—0 



N 



(4) / n ...^(r,( 



— ^ t 0-i 1 ...i N ijk T 
i,j,k=0 



i+j+k 



x sm* + *$ cos k 9 cos* 6 sinP 



Having specified the tensor components /^...i^ with re- 
spect to the Cartesian frame we derive the components 
related to the local orthonormal frame of spherical coordi- 
nates by a non-singular coordinate transformation. 

In order to infer the analytical properties of a scalar func- 
tion it is useful to rearrange the sum in (Q) . We first collect 
all the terms referring to cosmcj) and sinm^ respectively. 
We write 



(5) /(r,M) 



M 



m—0 



(a m (r,9) cos m4> + b m (r, 9) sinmcj 



where a m (r,9) and b m (r,9) behave identically in the fur- 
ther procedure. We opt for cos 19 and sinW as basis func- 
tions in 9 which allows the application of FFT-techniques 
for this transformation. An immediate conclusion from 
dJ) and the case i+j even is that the coefficients a2 m (r, 9) 
and 62m (r, 9) have to be expanded in terms of cos 16 while 
from the odd case that the expansion of a2 m +i(r,0) and 
&2m+i(?", 9) has to be done on the set sin/0. We therefore 
specify 



(G) 



(7) 



a,2m(r,9) = a li2m {r) cosW, 



a 2m+ i(r, 9) 



1=0 



(r) sin 19. 



Note that due to the well defined parity of aim, and a^m+i 
these coefficients — a priori only defined for < 9 < ir - 
can be continued analytically to periodic functions of 9 on 
the interval [— tt, it]. 

In the same manner as before we find from (^) that the 
polynomials di(r) are symmetric with respect to the inver- 
sion r — ► — r for I even and antisymmetric for I odd. As 
basis function set in r we decide for Chebyshev polynomi- 
als due to their superior properties in finite approximation 
schemes of non-periodic functions and the availability of 
fast Chebyshev transforms. Simplifications of the expan- 
sion scheme in the presence of additional symmetries are 
precised in the following paragraphs. 

Any regular function admits an expansion of this kind, 
but regularity furthermore implies additional constraints 
on the different coefficients. Having set up regular initial 
data according to (||), regularity of the involved quantities 
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is maintained during a calculation by the application of 
regular operators — we here ignore the influence of nu- 
merical effects due to aliasing or roundoff errors. 

Let us consider the covariant derivative of a vector in sphe- 
rical coordinates. For the scalar potential U(r, 8, <fr) = 
r sin 9 cos 6 and 



x sin 2Z 6> cos 2m 6> cos 2l ~ k d) sin k c 



(8) 

we have 



1 



U l% = {d r U, - d e U, 



1 



r sin i 



d+U) 



(9) (U r ,Ug,U<p) = (sin 9 cos (f>, cos 9 cos (f>, — sin < 
The covariant derivative Uo\a which reads 



(10) 



Ue\<f) — 



r sin0 

can be rewritten to yield 
1 



dd, Ug - 



r tan# 



(11) 



Ug 



r sin f 



{ds Ug - cosOUa). 



A numerical evaluation according to (|8j) and (|ll| ) reveals 
a perfectly regular behaviour. 

2.2 Supersymmetric case 

Additional spatial symmetries of physical systems involve 
continuous symmetry operations like rotation about a dis- 
tinct axis with an associated Killing vector field or discrete 
transformations like inversion at the equatorial plane 2 = 0, 
leading to distinct parity properties of the different tensor 
components. 

We define supersymmetry Jl9[ by the following behaviour 
of a scalar function: / is invariant with respect to inver- 
sion at the z— axis, hence f{—x,—y 1 z) — f{x,y,z), while / 
is [anti-] symmetric with respect to reflection at the equa- 
torial plane z = 0, hence f(x, y, —z) — ±/(x, y, z). Conse- 
quently / can be expanded into a sum of the type 



(12) f+(x,y,z) 



L 21 M 



2^ 2^ 1^ c khnX 21 k y k z 2m even case 

(=0 fc=0 m=0 



L 21 M 

(13) f-(x, V,z)=J2J2J2 cki m x 2l ~ k y k z 2m+1 odd case 

;=o fc=0 m=0 

which defines subsets of the general scalar functions intro- 
duced in Sec. 2.1. A write-up in terms of (r, 9, (f>) gives 



(14) 



L 21 M 

EEE 

1=0 k=0 m=0 



Cklm r 



2(l+m) 



(15) /-(r,< 



L 21 M 

EEE 

1=0 k=0 m=0 



Cklm r 



2(;+m) + l 



x sin 2i 6> cos 2m+1 6> cos 2l - k (h shA 



which can be modified, replacing sin 2 # by (1 — cos 2 



L 21 M 



r 2(l+m) 



(16) / + (r,M) = EEE 

Cklm ' 

1=0 k=0 m=0 

x(l-cos 2 6') i cos 2 

L 21 M 

(17) />,M) = EEE c ^ 2( ' +m)+1 

1=0 k=0 m=0 

x(l-cos 2 0)' cos 2m+1 <9 cos 2l ~ k (f> sinV 

According to Sec. |2.1| we conclude for the decomposition 
of a supersymmetric function: 

1. / is 7T— periodic in the angular variable 4>, 

2 • a2m(f,^) has to be expanded into a sum of cos 18 where 
I is even in the symmetric case and odd in the anti- 
symmetric one, 

3. ai(r) is an even polynomial in r for / symmetric and 
an odd polynomial for / antisymmetric. 

Physical problems which imply the use of supersymmet- 
ric functions allow to restrict the computational domain 
to [0,7r] in <fi and to [0,7r/2] in 9 which leads to an overall 
reduction of the effective grid size by a factor four. 
For the components of a vector field associated with the 
case of even supersymmetry we conclude: U r can be ex- 
panded into a sum of cos 21 9, Ug into a sum of sin 21 9 
and Us into a sum of sin(2Z + l)# which means that Ug 
undergoes a change of sign by reflection at the equatorial 
plane while U r and U& remain unchanged. The expansion 
of the radial part is done in terms of odd powers of r for 
all components. 

2.3 Axisymmetric case 

Axisymmetry restricts the set of scalar fu ncti ons under 
consideration to functions according to Sec. 2.1 where the 
only remaining term in (^) is ao (r, 9) . Thus / can be writ- 
ten as 

L 

(18) /(r,fl)=5^Si l0 (r)cosW 

1=0 

where az,o( r ) is an even function in r for / even and an odd 
function for I odd. 
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For the components U r and Ug of a vector field compatible 
with the assumption of axisymmetry we conclude: U r has 
to be expanded in terms of cos 19 where the radial part is 
even for I odd and odd for I even. Ug has to be expanded 
in terms of sin 19 where we also have a parity change in r. 
The properties of the lacking component which we need 
to handle in the models of rotating neutron stars can be 
derived from the Killing equation linked to axisymmetry. 
A short examination reveals that has to be expanded 
in terms of sml9 with a parity change in r — it behaves 
identically as Ug. 

2.4 Spherically symmetric case 

For sake of completeness we add the spherically symmetric 
case. As a further restriction of the axisymmetric case we 
keep from (|l§|) only ao,o(?")- f( r ) therefore reads 



K 



(19) 



/(*■) = E 



a 2 kr 



..2 k 



k=0 



which is an ordinary even polynomial in r while U r is rep- 
resented by an odd polynomial in r. 

3 Spherically symmetric neutron 
star collapse 

3.1 Basic equations 

The investigation of neutron star equilibrium configura- 
tions in spherical symmetry was the first problem in the 
fully general relativistic regime being solved by us by 
means of the (3+1 ) -formalism of general relativity and a 
spectral method 0, || ||] . A favourable choice of the line 
element ds 2 = g a p dx a dx@ in the case of spherical symme- 
try is given by RGPS {Radial Gauge-Polar Slicing) coor- 
dinates |20f| and reads 



(20) ds 2 = -N 2 dt 2 + A 2 dr 2 



{d9 2 



sin 



2 ). 



We stress the particular nature of this problem where the 
field variable A is not really a dynamical quantity. In fact it 
is uniquely determined at any moment as well as the lapse 
function N by the matter fields which have to be evolved by 
means of the hydrodynamical equations. Since the solution 
outside the star is known in advance to coincide with the 
static Schwarzschild solution of a point mass of the same 
size according to the Birkhoff theorem, we benefit from a 
double simplification. First the whole time evolution is yet 
determined by propagating the hydrodynamical variables 
and further the calculation can be restricted to the stellar 



interior. The interior solution for the gravitational field 
has then to be matched to the analytical exterior one. We 
further note that due to the static character of the exterior 
solution no gravitational waves — which otherwise would 
be an observable of most importance — are emitted. 
Concerning the hydrodynamical part we employ a set of 
particular variables which lead to equations ressembling 
very closely their Newtonian counterparts including a gen- 
eral relativistic generalization of the classical Euler equa- 
tion. We note the privileged role of the Eulerian or local 
rest observer Oq in this formulation. The hydrodynamical 
equations read 



(21) 



1 



d t £ + —d r (r 2 {£+p)V r ) =0, 



(22) dtU + V r d r U 



1 



N 



d r p + Ud t p 



£+P \A 
AN fm 



(23) 
(24) 



d t D+-^d r {r 2 DV r ) = 
dtSs + V r d r SB = 



where we have introduced the energy density £ and the 
fluid velocity U measured by Oq, the coordinate baryon 
density D and the entropy per baryon sb while V r de- 
notes the fluid coordinate velocity. The Lorentz factor r 
is defined as r= (1 — U 2 )~z . In addition one has to solve 



(25) 



d r m{r, t) = Airr 2 £{r, t), 



(26) <9 r $(r, t) = A 2 



m{r, t) 



+ 4?rr[p+ {£+p) U 2 



where A(r,t) is related to m(r,t) through 

2m(r, t) 



(27) 



A{r,t) = 1 



The neutron star matter is modeled as a perfect fluid, 
adopting a realistic dense matter equation of state. 

3.2 Numerical method 

The initial value model for the dynamical calculations 
is provided by solving the Tolman-Oppenheimer-Volkoff 
equations describing a spherically symmetric static star 
where each model is determined by the central value of 
the pseudoenthalpy H c . A first solution is obtained by in- 
tegration of this system of ordinary differential equations 
while an overall numerical accuracy of the order of 10" 14 , 
adapted to the subsequent use of a spectral method, is 
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2 4 



t [ms] 

Figure 1: Relative variation of the stellar radius with 
elapsing time for a stable equlibrium model with < 
n% h and \AM/M max \ = -5xlCT 5 . The hydrodynamical 
timescale is indicated on the lower left. 




Figure 2: Profiles of the lapse function N for t ranging from 
to 7.296 ms. The location of the Schwarzschild radius R§ 
is indicated by a vertical bar. Direction of increasing time 
is downward. 



achieved by iteration of the approximate solution. 
The computational domain is identical with the stellar in- 
terior during the whole calculation, thanks to a comoving 
grid whose outer boundary coincides with the star surface. 
This maintains a constant spatial resolution during the col- 
lapse and a fine sampling of the steepening gradients near 
the star surface thanks to the accumulation of the Gaufi- 
Lobatto points in this region. It furthermore minimizes 
the advective terms, hence improving the numerical accu- 
racy. All quantities are expanded in terms of Chebyshev 
polynomials in r mapping the interval [0, R*(t)] onto [0, 1] 
taking into account their analytical properties according to 



Sec. 2.4. The 2 nd order semi-implicit time integration en- 



sures the stability of the code which allows us to perform 
simulations of a duration of many dynamical timescales. 
This ability is very important when studying the effects 
of perturbations on equilibrium configurations. This task 
is also favoured by the fact that we integrate the original 
system of equations without any artificial viscosity to sta- 
bilize the code while in addition the intrinsic viscosity of 
spectral methods is negligible. The ingoing characteristic 
of the hydrodynamical system at r = i?* gives rise to one 
boundary condition which is chosen to fix the baryon den- 
sity at the star boundary. It is imposed by means of a 
r-Lanczos procedure 21 on the system as a whole which 
is the well posed mathematical procedure. 



3.3 Results 

Among the various tests we have imposed on our code one 
describes the collapse of a homogeneous dust sphere whose 
solution was given by Oppenheimer and Snyder . From 
the beginning of the collapse until the moment where the 
whole configuration is highly relativistic and practically 
frozen we observed the different variables to reproduce the 
analytical values within an error of better than 1CT 5 @, 
while the errors related to previous studies based on finite 
difference methods are of the order of 10 -2 [ p3| . 
Neutron star models near the maximum mass configura- 
tion — M max = 1.924M Q and R = 10.678 km for the em- 
ployed EOS — are interesting with respect to their sta- 
bility against radial perturbations. Configurations with a 
central baryon density approaching the critical one exihibit 
an oscillatory behaviour which is dominated by the funda- 
mental mode of oscillation of increasing period length. It 
represents a uniform growing and shrinking of the entire 
star which is modulated by less important harmonics of 
higher order. Fig. |l| shows this temporal variation for a 
stable configuration. We stress that these oscillations are 
entirely driven by roundoff and discretization errors of a 
total order of 10~ 10 — no external force has been applied 
to trigger this variability. Increasing the central baryon 
density beyond the critical value one enters the branch of 
unstable configurations. This time the fundamental mode 
starts a contraction of the star which results in an unlim- 
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Figure 3: Relative variation of the total baryon number 
during the collapse. 

ited collapse. Though the actual coordinate choice is not 
capable to properly capture the formation of an appar- 
ent horizon which would clearly reveal the formation of a 
black hole the evolution of metric potentials and matter 
variables gives a distinct indication of this event. Take for 
instance Fig. || which shows the time development of the 
lapse function N, measuring the elapsed proper time of 
the local Eulerian observer Oq. Inside the Schwarzschild 
radius it tends toward zero for increasing coordinate time 
t. This behaviour — called the 'lapse of the lapse' - 
is related to the singularity avoidance property of the cho- 
sen coordinate gauge which keeps the spatial hypersurfaces 
from propagating into a forming singularity. The dynam- 
ical timescale of the collapse is the time elapsed from the 
beginning of the collapse until the moment where the evo- 
lution appears to be frozen to a distant observer and has 
the value t — 7.4 ms for the considered configuration. Fig. || 
illustrates the relative error committed on the total baryon 
number which is a conserved global quantity. In the early 
phase it is conserved with a relative accuracy of 2 x 10 -8 
and with 4xl0~ 6 during the violent stages of the collapse. 
The deviation increases up to 5 x 10~ 5 in the final phase 
where sharp gradients form near the horizon. The proper 
working of the code in the perturbative regime has been 
recently confirmed by direct comparison with a linear adi- 
abatic code Q . The calculations have shown a very good 
agreement of the frequencies of the fundamental modes for 
the two opposite approaches. The original version of our 
code has been extended later on to include the processes 
of neutrino production and transport during neutron star 



collapse in order to compute the observable neutrino emis- 
sion for a distant observer [ fli"| . A multidomain extension 
of this code is in use in order to simulate type II super- 
novae in spherical symmetry. It is particularly well suited 
to properly capture the high contrast up to ~ 10 6 of the 
mean densities in the individual subdomains which appears 
during the collapse. Four to five zones are typically needed 
to cover the dense central core forming the new born neu- 
tron star and the outer layers of much lower density with 
the desired resolution. In contrast with the original code 
boundary conditions are imposed in this case by means of 
a modified r-Lanczos scheme applied in coefficient space. 

4 Axisymmetric rotating relativi- 
stic bodies 

4.1 Basic equations 

In order to study rapidly rotating neutron stars we have 
made the assumptions of spacetime being axisymmetric, 
stationary and asymptotically flat. We have further sup- 
posed spacetime to be circular, thus the absence of merid- 
ional currents in the sources of the gravitational field. In 
this case spacetime can be described by MSQI (Maximal 
Slicing- Quasi Isotropic) coordinates jl3j which have a line 
element ds 2 = g a p dx a dx@ of the form: 

(28) ds 2 = ~N 2 dt 2 + A 4 B- 2 (dr 2 + r 2 d6 2 ) 

+A 4 B 2 r 2 sin 2 (9 (# - N^dt) 2 . 

Spacetime is hence fully determined by the four metric po- 
tentials N , N't', A and B. MSQI-coordinates are global 
coordinates and lead to elliptic operators which admit a 
consistent treatment of boundary conditions and an effi- 
cient solution by spectral methods. The matter fields are 
chosen to model a perfect fluid where first a polytropic, 
hence analytical, equation of state was used. The assump- 
tion of a perfect fluid reduces the equations of motion to 
an algebraic equation for the heat function H. From the 
Einstein equations one derives four elliptic equations for 
the variables v = h\N , N^ and 

(29) G(r,6) = N(r,9)A 2 {r,9)B(r,9), 

(30) C{r, 9) = v{r, 9) + 2a(r, 9) - (3{r, 9) 
whose final form reads 

A 4 

(31) A 3 v = — 7 [4TT(E+S l t ) + 2(k 2 1 +k 2 )}-dvd(is+2a+p), 

(32) A 3 A> = -167r^^^-rsin6>GW d(6cH-3/?-i/), 

B 4 r sin 9 
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(33) 



(34) 



N A 6 

A 2 G = 8Tr——rsine(S r r + S e e ), 
H 

^S^^ + i{kl + k 2 2 )]-{du) 2 



A 2 C 



^1 

B' 2 



where we have introduced a — hi A, [3 — h\B and 

(35) G(r,0) = rsin0G(r,0), 

(36) N*(r,0) = r sin 0N^{r, 0) 
as well as the abridged notation 



(37) 



da 3(3 — d r a d r (3 — - > 



9ade(3. 



Further employed quantities are the total energy density 
E, the stress tensor S t j , the momentum density Ji and k\ , 
fc 2 which are related to the extrinsic curvature tensor Kij . 
A2, A3 and A3 denote scalar Laplacians in two and three 
dimensions and a vector Laplacian in three dimensions re- 
spectively. 

Existence and uniqueness of the solution of these ellip- 
tic equations are ensured for physically relevant cases 

nil Hi- 

Note that ( |33| ) and (|34|) can be continued analytically to 
yield genuine 2D-Poisson equations in the entire (r,6)— 
plane. One infers immediately that £ exhibits a logarith- 
mic divergence for r — > 00 unless the total integral over 
the source of ( |34| ) vanishes identically whereas we require 
(\r=oo = 0. This 2D-virial theorem of general relativity 
(GRV2) [|l6] is in some sense related to the classical New- 
tonian virial theorem and furnishes a consistency condition 
for any solution of the Einstein equations which is compat- 
ible with our basic assumptions. It has to be taken into 
account during the calculation and provides a strong con- 
sistency check of the numerical solution. 

4.2 Numerical method 



The mathematical problem involves (31 )-(|34|) and an al- 
gebraic first integral equation for the matter fields. Our 
numerical solutions are exact in the sense that the gov- 
erning equations are derived from the full theory of gen- 
eral relativity without any analytical approximation while 
the numerical code solves these equations in all space ex- 
tending the numerical integration to spatial infinity which 
allows to impose the exact boundary conditions of asymp- 
totical flatness on the gravitational fields as well as the 
proper calculation of the source terms of (|3l|) - (|34|) which 
fill all space. 

This is accomplished by the use of two grids, where the first 
one covers the stellar interior using the radial variable r in 



the interval [0, R], while the outer space is compactificd 
thanks to the variable transform u = r~ l and in this way 
mapping [R, 00] onto the finite interval [i? _1 ,0]. Whil e in 
the 0— variable a Fourier expansion according to Sec. |2.l| 
is used, the radial part is expanded in terms of Chebyshev 
polynomials. The inner zone is mapped onto half the def- 
inition interval [0, 1] which allows to take into account the 
parity properties of regular functions with respect to the 
origin according to Sec. 2.3. In the compactified zone the 
expansion has been done in the usual manner on the whole 
definition interval [— 1, 1] of Chebyshev polynomials. 
The effective scheme which is based on a relaxation method 
works as follows. We consider rigidly rotating neutron 
stars with a polytropic equation of state. A particular 
configuration is hence determined by fixing the value H c 
of the heat function at the centre of the star and its angu- 
lar velocity f2. We start from very crude initial conditions 
where all the metric quantities are set to their flat space 
values (a, j3, v, £ and — G—l) and the matter dis- 
tribution is determined by a first approximate guess. 



While the GRV2 identity related to (34) holds for an exact 



solution of the Einstein equations we have to enforce this 
consistency relation at each iteration step in order to avoid 
a logarithmic divergence of the approximate one. This can 
be accomplished by modifying (b3) according to 



(38) A 2 £ = at + a 



A 2 £ = ^ 



C 



(39) al = 8^S%, a« = ^{3(kl + k 2 2 )]-(dvf. 

At each iteration step A is chosen in such a way that the 
total source integral vanishes. The final solution has to 
satisfy d34| ) exactly which is equivalent to A = 1 . The de- 
viation of A from unity during the iteration measures the 
violation of self-consistency of the approximate solution 
and can be used to monitor convergence. 

The sources of (H|) -(H) exhibit some terms involving 
simple operators like r, smO, d r , etc. which are accurately 
computed in coefficient space before evaluating the entire 
expressions in configuration space. Expansion of the total 
sources in terms of the angular eigenfunctions of the differ- 
ent Laplacians (P; O (cos0), P^(sin0) and (cos 10, sin 10) for 
A 3 , A 3 , and A 2 respectively) leads to a system of ODEs in 
the radial variable. The unique global solution is obtained 
by appropriate linear combinations of the corresponding 
particular and homogeneous solutions in order to match 
the piecewise solutions at the grid interface and to satisfy 
the boundary conditions at r = oo. The new values of the 
gravitational field variables are then used to update the 
matter distribution by means of the first integral equation 
and the iteration can go on. 
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Figure 4: Comparison between numerical and analytical 
solution for the Schwarzschild interior and exterior incom- 
pressible solution. The location of the star surface is indi- 
cated by an asterisk. The plotted quantities are the relative 
error in the pressure with respect to the central pressure 
inside the star and the absolute error in N (solid line) and 
A 2 (dashed line) outside the star. 

4.3 Tests 

We have subjected our code to two different kinds of tests 
which ensure the reliability of the numerical results. 
External tests consist in the comparison with previous so- 
lutions, either analytical or numerical ones. Such a test of 
the code has been performed for an analytically known 
Schwarzschild type solution of a non-rotating homoge- 
neous sphere with the corresponding numerical one. Rel- 
ative errors committed on global quantities such as total 
gravitational mass and circumferential radius are of the or- 
der of 1CP 14 . This accuracy holds also for local quantities 
as shown in Fig. |^ for the pressure p inside the star and 
the metric coefficients outside the star — none of the errors 
exceeding 10~ 14 . A recent project of systematic calibra- 
tion and comparison of the numerical results of different 
groups working in this field has yielded an agreement of 
characteristic quantities of realistic neutron star models at 
a level of about 1CP 3 . 

Internal tests represent the second important class of tests 
and are derived from some relations of global or local char- 
acter which are related to the Einstein equations but not 
automatically enforced during the calculation. These tests 
are very powerful, since they do not only verify the proper 
working of the numerical scheme for some simple — usually 




-10 10 

x [ km ] 



Figure 5: Level contour E(r, 9) in the case of a polytropic 
EOS with 7 = 2 for fl = n K . 

degenerate — test problem, but apply to any calculation 
and supply an intrinsic estimate of the numerical error in- 
volved. In the following neutron star matter was modeled 
by a pefect fluid with a 7 = 2 polytropic EOS. The angular 
velocity has been varied between f2 = for the static case 
and Q = f2x for the maximum rotating case where for a fur- 
ther increase of mass shedding along the equator occurs. 
Fig. H shows the flattened shape of a neutron star rotating 
at f^K- A simple test makes use of the known analytical 
expansion of B according to B = l+r 2 sin Of. It showed 
B to coincide with the analytical value 1 on the polar axis 
within 1CP 6 in the most unfavourable case of maximal an- 
gular velocity. The principal test is provided by the GRV2 
identity. An examination of the Schwarzschild type solu- 
tion has revealed that |1 — A| is very closely related to the 
global errors derived from the numerical solution for vari- 
able N r . This observation, though obtained for the static 
case, is supposed to hold in the rotating case as well. |1— A| 
can thus be considered as an estimator of the global nu- 
merical accuracy. Figs. ^], [7] show the convergence of 1— A| 
during the iteration for S7 = and f2 = f^K ■ While in the 
first case the exponential decay of |1 — A| continues until 
roundoff errors of 1CP 14 mark a lower limit, the total er- 
ror in the rotating case is about 10~ 6 . This difference is 
due to the deviation of the flattened stellar shape from the 
spherical numerical grid which leads to a discontinuity in 
the derivatives across the stellar surface located inside the 
inner zone, where the diverse quantities hence are no more 
analytical functions. The attainable accuracy in depen- 
dence of the number of grid points is illustrated in Fig. ||. 
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Figure 6: Convergence of the error indicator |1 — A| during 
the iteration process for f2 = 0. 



Figure 8: Internal error indicator as a function of N r in the 
compressible spherically symmetric case (polytropic EOS). 
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Figure 7: Convergence of the error indicator |1— A| during 
the iteration process for fl — £1k- The spike at N = 10 is 
due to switching on the rotation. 



The exponential decrease, usually called evanescent error 
and characteristic for spectral methods, is clearly visible. 
This property is lost in the rotating case where we observe 
a power law decay of the committed error oc N~ 4 - 5 as found 
from Fig. |j. This inconvenience will be overcome by the 
implementation of an adaptive ellipsoidal grid which aligns 
the domain boundary to the star surface. We finally con- 
clude that we have computed neutron star models with an 



Figure 9: Same as Fig. || but for maximal angular veloc- 
ity 11k- Pay attention to the log — log scale; the spectral 
properties are lost here and one observes a power law con- 
vergence. 

analytic EOS achieving a precision of 1CP 14 in the static 
case to some 10~ 6 in the maximum rotation case which 
has to be compared with previous results based on finite 
difference methods of the order of 1CT 2 ||, || §|, HJ- 
It is further interesting to note that for typical values of 
N r = 32 and Ng = 16 one iteration is performed in 480 ms 
on a VAX 4500. For an average number of 50 iterations per 
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model the whole calculation is finished in about 24 s. The 
code efficiency has enabled us to carry out extensive stud- 
ies of neutron star samples under employment of numerous 
realistic EOS of neutron star matter [Q. 

4.4 Rotating neutron stars with magnetic 
field 

A further step toward a realistic description of rapidly ro- 
tating neutron stars has been recently achieved by the fully 
self-consistent inclusion of magnetic fields into our models 
These calculations represent the first numerical so- 
lutions of the coupled 2D-Einstein-Maxwell equations for 
rotating neutron stars which are hence fully relativistic, 
taking into account any kind of interaction of the electro- 
magnetic field with the star and the gravitational field. 
To complete the physical specification of our neutron star 



Relative error on electric field 



Relative error on magnetic field 



models as described in Sec. 4.2 we assume a perfect conduc- 
tor behaviour of neutron star matter — the star interior is 
thus free of electric fields — and add the electromagnetic 
field variables A t , A$, the current variables j t and j$ and a 
structure function / which determines the current distribu- 
tion inside the star. The additional free parameter which 
fixes a unique neutron star configuration is the total elec- 
tric charge Q. A derived global quantity is the magnetic 
dipole moment M. which characterizes the magnetic prop- 
erties of the neutron star. The Maxwell equations lead to 
a set of coupled elliptic partial differential equations which 
exhibits a similar structure as j3l|)-(|3^). They involve a 
scalar Poisson equation for A t and a vector Poisson equa- 
tion for Aa, which read 



(40) A 3 A t 



A 4 



A i B 1 
N 2 



A V sin 2 9 x dA t dN^ 



- 1 



(r sin 



(41) A 3 i* 



A i B 2 

~{dA t + 2N' t >dA ! f ) ) d(2a+[3-v) 

f 1 \ 

-2 drAj, + — — - deA* , 

r \ rt&nV J 

-fi A s (/-jVV)rsin0 
A^B 2 



N 2 
1 



r sin 9 dN 4 '(dA t +^8 A<f,) 
dA (j) d{2a+(3-v) 



where we define 
(42) 



A,i 
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Figure 10: Comparison between the exact and the numer- 
ical solution in the case of a rotating conducting sphere 
with a magnetic point dipole at its centre. The asterisks 
indicate the sphere's surface. 



In this extended formulation the determination of A t and 
A0 precedes the solution of the former set of equations. 
A<p is a smooth function in all space — we recall that 
fi r ~ 1 for neutron star matter — and can therefore be 
uniquely solved by imposing the asymptotic boundary con- 
dition A$ — > for r — > oo. The treatment of A t is slightly 
more complicated due to its non- smoothness across the 
stellar surface — a behaviour caused by the surface charges 
which are characteristic for ideal conductors. Since A t is 
linked to A^ by a linear relation in the star interior, the 
exterior solution has to be matched to the latter one, re- 
specting the condition At = at infinity. We note that, 
once more, this is a task which is easily performed thanks 
to the use of a spectral method. After solution of the whole 
system we can update the electromagnetic contributions to 
the stress-energy tensor and proceed with the solution of 



the Einstein equations according to Sec. 4.2 



r sin 6 



Test calculations of the electromagnetic part have been 
performed for some simple cases, where analytical solutions 
are known, among these one involving a rotating magnetic 
dipole where an infinitesimal current loop at the origin is 
surrounded by some rotating perfectly conducting sphere. 
This testbed calculation mimics the configuration of a ro- 
tating neutron star with surface charges. Fig. |l0| shows 
the relative error committed on the electric and the mag- 
netic field for a conducting sphere of 1 m radius rotating 
at — 3000 s _1 with a central current loop of jo = 10 11 
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Figure 11: Comparison between Ferraro's analytical so- 
lution and the numerical one in the case of a Newtonian 
incompressible fluid endowed with a magnetic field corre- 
sponding to a constant current function /(&) = fa- The 
plotted quantity is the relative difference between the two 
values of as a function of the radial coordinate r for 
three values of 6. Asterisks indicate the star surface. 



Am -2 . Apart from the origin where the model current dif- 
fers from the theoretical 5-distribution the relative error is 
very small: about 10 -5 at half the radius and reaching its 
minimal value of 10 -9 outside the sphere. An analytical 
solution for a Newtonian incompressible fluid J3l| endowed 
with a particular current distribution under the — simpli- 
fying — assumption of spherical symmetry was adopted as 
a more sophisticated test. Also in this case the agreement 
was quite good with a relative error of better than 10~ 3 , 
shown in Fig. |ll|. The deterioration with respect to the 
dipolc problem is due to some simplifying assumptions of 
the analytical model. 

The accuracy of solutions of the complete Einstein- 



Maxwell equations was estimated as in Sec. 4.2 by use of 
the virial identity GRV2 JlTj as well as by a more general 
three dimensional integral identity valid for any stationary 
and asymptotically flat spacetime which we call GRV3 [p| . 
It is the general relativistic generalization of the classical 
Newtonian virial theorem. The actual values of 1 — A 
showed about 10 -5 for analytical EOS and some 10 -4 for 
the tabulated ones which agree with those of calculations 
without magnetic field. Throughout all the calculations 
we chose a grid resolution of 41 points in r and 21 points 
in 6. 
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Figure 12: Fluid proper density isocontours in the (r,9)— 
plane for the M = 4.06M© maximum mass static magne- 
tized star built upon a polytropic EOS for 7 = 2. The thick 
line indicates the star surface. 




km 



Figure 13: Magnetic field lines in the (r,8) -plane for the 
maximum mass configuration corresponding to Fig. |l2|. 
The thick line indicates the star surface. The magnetic 
field amplitude amounts to Be = 9 x 10 4 GT at the star's 
centre. 



In the following we studied configurations of static neutron 
stars endowed with a magnetic field. The absence of kine- 
matical effects allows an unambiguous interpretation of the 
effects of the magnetic field. In the static case the electric 
charge vanishes identically, leaving alone a magnetic field. 
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Figure 14: a. Magnetic field lines in the (r, 9) -plane for the 
configuration specified in the text. The thick line indicates 
the star surface, b. Electric isopotential lines A t =const. 



lows us to adopt elliptic equations for the gauge variables 
N and N l which have a direct geometrical signification. 
The maximal slicing-minimal distortion gauge ^J, Q is 
highly singularity avoiding and neatly captures the propa- 
gation of gravitational waves in the far field zone. This is 
a very important feature, since the study of gravitational 
wave emission is the principal goal of this investigation. 
The governing equations then have the following form: 

(43) d t hij +Ni\j + Nj | i + 2NKij = , 

(44) d t K l J +N l K l h i-N l .iK l j + N\ 1 K l i 
+N\\ j -N[R l j +An{S-E) ^ i -87rS" i ] = 0, 



Since the stress-energy tensor is not isotropic, we observe 
a deformation of the star already in the static case. Fig. |lj 
shows a maximum field configuration and Fig. [l^ the cor- 
responding distribution of the magnetic field. In the ro- 
tating case the magnetic field is accompanied by an addi- 
tional electric field. The both field distributions for a Pol2 
M = 3.37M Q model at Q — 3 x 10 3 rad s _1 and a constant 
current function are given in Fig. [l4|. Note that the field 
lines of A t and A<p coincide in the star interior due to the 
perfect conductor assumption. The non-smoothness of A t 
across the star boundary properly reflects the discontinu- 
ity of the electric field due to the existing surface charges. 



5 3D— Gravitational neutron star 
collapse 

5.1 Basic equations 

The investigation of the 3D-gravitational collapse of rotat- 
ing neutron stars requires the solution of the general time- 
dependent field equations of general relativity. We stress 
that already in the Newtonian case a fully three dimen- 
sional simulation of stellar collapse is a highly demanding 
and non-trivial problem. Indeed up to this day there ex- 
ists only one corresponding investigation which aimed at 
the study of gravitational wave emission associated with 
type II supernovae J3^| . While in the presence of symme- 
tries like stationarity or axisymmetry the field equations 
are greatly simplified when evaluated explicitly for some 
appropriate coordinates the situation is quite contrary for 
a three dimensional dynamical problem. It is therefore 
favourable to solve the Einstein equations in their original 
coordinate independent (covariant) form. This choice al- 



(45) R-ICj K j i-16irE = 0, 

(46) 871-^ = 0, 

(47) d t E+N l E tl ~NK l 3 S 1 l +N~ 1 (N 2 .r) ll = 0, 

(48) d t J % +N l J^ l +N\ t Ji + {E5h + Sh)N\ J +N S 3 Aj =Q, 



(49) 



N\ l ll -N[K l ] K j i + 4Tr(E+S)}=0, 



(50) JST» b - + - {N j u Y + WjN 3 + 2 K ij N {j + WitNJ 1 = 0. 
o 

This system of coupled partial differential equations is 
characterized by the following properties: 
The system includes a time first order hyperbolic system 
(43) and (Q) of the spatial metric tensor hij and its con- 



jugate momentum variable Kij which reduces to a wave 
equation for hij in the far field zone. The evolution of the 
matter fields E and Ji is governed by the parabolic system 
(47) and (ff8|). These equations constitute the dynamical 
part of our problem. N and N z are subjected to the Pois- 
son type equations ([l9|) and (|5(]) where the matter fields 
act as source terms. Further involved quantities are the 
stress tensor Sij and the Ricci tensor Rij where S = S l i 
and R = R l i. The additional constraint equations (|45|) and 



(46) which establish some consistency relations between 



gravitational and matter fields are satisfied identically for 
the exact solution. They may be used to reduce the num- 
ber of the dynamical variables, an approach that would 
result in a constrained evolution scheme. Indeed this pro- 
cedure guarantees that the numerical solution represents 
at any moment some solution of the Einstein equations 
but not necessarily the correct one. We prefer an uncon- 
strained scheme where the constraint equations can serve 
to estimate the involved numerical errors. 
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5.2 Numerical method 



In the spirit of the analytical approach of Sec. 5.1 we have 



opted for a numerical scheme which is based on a one-to- 
one adaptation of common tensor calculus, where elemen- 
tary operations like contraction and covariant derivation 
are performed by specific subroutines acting on entire ten- 
sor quantities. We further introduce a flat background 
metric which enables us to separate the contributions re- 
lated to the curvature of space and to carry out the nu- 
merical operations with respect to flat space spherical co- 
ordinates and the corresponding metric tensor /y p3|] . For 
a given three-metric hij we then have the following rela- 
tions: 

We first introduce a tensor A z jk defined as 

(51) A l jk = - h ll (hji\\k + hki\\j - hjk\\i) 

where || denotes covariant differentiation with respect to 
flat space spherical coordinates. The covariant derivative 
U l \j of a vector U l in curved three-space can then be re- 
expressed as 

(52) U i ]j = U i h +A i jk U k 

where the generalization to tensors of higher order is obvi- 
ous. We can further rewrite the Ricci tensor Rij in terms 
of the A l jk ■ 

(53) — A 1 — ^ + A / yA m / m — A' i m A m ji. 

The effective use of our scheme is clarified by inspection of 
the already familiar scalar Poisson equation. 



(54) iV 1 



* m ijN\\ 



= 5. 



Defining a new tensor field r as = f^+r^ we can isolate 
the (/y)-related covariant Laplacian and find the bimetric 
equivalent of the covariant scalar Poisson equation 



(55) N^ l{i = S+(f ij +r ij )A' 



In the case of a conformally flat metric — A 4 /y with 
fij being the usual metric tensor of flat space spherical 
coordinates we obtain r y = (A -4 — 1) and the following 
equation where Af represents the usual flat space scalar 
Laplacian (a denotes In A). 

(56) A f N = S -2A- i (d r ad r N + -^d e ad e N 



The linear contributions on the right can be eliminated, 
if we slightly change the definition of r to describe the 



deviation of the conformal metric h — r y~ 1 h from the flat 
space metric /, where 7 3 = | h 1 1 / | ~ 1 . So with h 1 ^ = fv+r^ 
we have r — and 



(57) A f N = A 4 S - 2 (d r a d r N 

1 



+ \ dea d$N 
Y-d^ad^N). 



Though this result has been derived for a conformally flat 
metric we conclude that for any three-space the separation 
of the conformal factor is a first improvement compared to 
the flat space approximation. 

Following our reasoning in Sec. |2.l| we restrict ourselves to 
the exclusive use of pseudophysical components of tensor 
quantities related to the standard local orthonormal frame 
of flat space spherical coordinates in our numerical scheme. 
The subroutines which calculate the covariant derivatives 
of tensors of order zero to three appearing in our equations 
show the typical relative errors of the order of 1CP 14 due 
to roundoff erors when using simple test functions. For a 
successive derivation of a scalar function down to a ten- 
sor of order four the relative error still does not exceed 
10 -12 . While the computation of the lapse equation can 
be reduced to the iterative solution of a scalar Poisson-like 
equation as demonstrated in Sec. |4.2| , there remain two 
other problems of higher demand. The first one concerns 
the solution of the general shift vector equation involving 
a linear vector operator which comprises a vector Lapla- 
cian and the gradient of a divergence applied to the shift 
vector N l . In Sec. 5.3 we present a decomposition scheme 



based on the Clebsch-Gordan theorem which leads to an 
equivalent system of three scalar Poisson equations that 
can be solved successively. The other one is related to the 
semi-implicit time integration of the evolution equations 
of the metric potentials. Here Rij includes a tensor Lapla- 
cian Ahij which has to be treated implicitly. An equivalent 
approach as in the vector case is in work. 

5.3 Vector Poisson equation 

The numerical inversion of a 3D-vector Poisson equation is 
necessary to solve the general shift vector equation which 
is an equation of the following type 



(58) 



AV + aV(V-V) = S 



where a is a constant. In order to facilitate the solution 
of ( |58| ) we decompose V and S into its divergence-free 
and its irrotational part. We introduce vector fields V 
and S which we suppose to be divergence-free as well as 
two scalar potentials W and <f>. We thus have a unique 
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decomposition 



where A f denotes the ordinary scalar Laplacian 



(59) 
(60) 



V = V 

s = s 



after specifying the appropriate boundary conditions for $ 
and In a first step we solve the Poisson equation for <f> 
which we obtain by taking the divergence of (p0|). 



(61) 



A$ = VS. 



Rewriting ( 58 ) in terms of the new variables we infer the 
equivalent equation 



(62) 



Ay + V((l+a) A* - $) = S. 



Taking the divergence of ( |62| ) where we make use of the 
commutativity of differential operators in flat space, we 
find a Poisson equation linking A 'J' to the already known 
potential <£>. 

(63) (l+a)AA* = A$. 

The solution of ( |63| ) is a priori determined up to an addi- 
tional potential $h with A$h = 0. We fix $h to be zero 
according to the boundary conditions and obtain a Poisson 
equation for "P. 

(64) (l+a)A* = $. 



We turn now to (|62j). Taking into account fl64| ) we derive 
the final equation that governs the divergence-free fraction 
of V. 

(65) AF = S. 

We specify the components V% to be the physical compo- 
nents of the vector V related to the local orthonormal basis 
of spherical coordinates. We further drop the tildes on the 
vector quantities. The explicit write-up of (p5J) then reads 



(66) 



AfVr oV r ~ d e 



1 



tan 

d A V A = S, 



V e 



(69) A f =d 2 r 



1 



r 2 tan 



do 



r 2 sin 2 ( 



We further recall the representation of the covariant diver- 
gence V • S in spherical coordinates 



(70) V-S = d r Sr + - S r + - deSfl 

r r 



1 



r tan 6* 
1 



So 



rsm.0 

We adopt the following notation where we introduce two 
auxiliary scalar potentials U and W in order to obtain a set 
of decoupled equations which is equivalent to (|66|)-(|68|). 



(71) 
(72) 



Vg : 



dgU 



d e W 



r sin ( 
1 



rs\n0 



From V • V = we get 



(73) d r V r + - V r 

r 



■dlu 



r 2 tan 



d e U 



dlU = 0. 



Combining (|73|) with (pq ) leads us to a scalar Poisson equa- 
tion for V where we have defined V = rV r . 



(74) 



AfV = rS r 



Proceeding in the same manner with ( p7|) and (|68J) we 
derive the following linear combinations of the resulting 
equations with ([T3h. 



(75) sm0do(d 2 U-d r V r ) - <9 ( AW ~ - d r W 

— r sin So , 

(76) ^— d^(d 2 U-d r V r ) + do [ AW - - d r W 
sm0 \ r 

= rSch 



(67) 



A f V e 



9 ■ 2 

r z sm 



(68) A f V^ 



r 2 sin tan i 



r 2 sin 2 ( 



■ i9<*V0 — So, 



r 2 sin 



r 2 sin tan 



Ba V r 



— Sa 



Derivation of ([75D with respect to and of (|76| ) with respect 
to and adding the both equations allows us to recover a 
scalar Poisson equation for (d r U— V r ) which only involves 
the angular variables (0, (j>) where we have integrated over 
r and set the implied integration constant to zero to assure 
a vanishing behaviour of the source terms at r = oo. 



(77) 



AeA(d r U-V r ) = -r 2 S r . 



Here Ao<p denotes the angular fraction of Af multiplied by 

^2 
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Figure 15: for a rapidly rotating Kerr black hole at 
a/M — 0.99 where rn denotes the radius of the horizon. 
The different curves correspond to various values of 9 rang- 
ing from 0° to 90°. The asterisks indicate the subdomain 
boundaries. 



At this point we can determine U by means of V r which has 
already been fixed by ( f74| ) . In order to calculate the lacking 
potential W we take up (|76|). An ordinary integration 
over 9 where we require vanishing behaviour of the source 
terms at infinity as for (d r U— V r ) results in a scalar Poisson 
equation for W defined by W = rW, 



(78) rAW = 



1 



d^dlU-drVr 



cW. 



This scheme hence allows to compute the required poten- 
tials successively for a given source distribution starting 
from the ordinary Poisson equation for V, solving then the 
equation involving (d r U — V r ) and thus fixing U while in 
a last step W can be determined as a quantity depending 
on V r and U. Note that Sg does not appear in the source 
terms of the final equations. This is due to the constraint 
equation V • S = removing one degree of freedom. 



5.4 Tests 

The routines we have built based on this computation 
scheme and on our spectral method library enable us to 
solve a vector Poisson equation in a multidomain config- 
uration including an exterior compactified zone — if de- 
sired — which covers all space and thus allows to impose 



Figure 16: Absolute error in corresponding to Fig. pi 



proper boundary conditions for asymptotically flat space 
where for the time being we ristrict ourselves to the su- 
persymmetric case. Numerical tests have been performed 
on simple test functions and on problems where the ana- 
lytical solution was known in advance. For test functions 
we found once more the numerical error to be governed 
by the roundoff limit of the employed machine of the or- 
der of 10~ 14 . As an advanced test problem we solved the 
shift vector equation in the Kerr metric of a rotating black 
hole. The presented configuration corresponds to a rapidly 
rotating Kerr hole — a/M ~ 0.99 — close to maximum 
angular velocity where the relativistic effects involving the 
shift vector component are strongly pronounced. The 
entire space outside the black hole is covered by three zones 
where the outer compactified one extends to spatial infin- 
ity. The grid resolution is chosen to be N r = 33 in each 
zone and Ng = 9. Note the quite small number of nodes 
in 9. Thanks to taking into account the even symmetry 
of the problem with respect to reflection at the equatorial 



plane according to Sec. 2.2 this corresponds to an effective 
value of 17 in the whole interval [0, ir}. Fig. [l5| shows the 
course of in the vicinity of the black hole while Fig. [l6] 
illustrates the absolute error committed on the shift vec- 
tor component . The numerical error nowhere exceeds 
10~ 10 . For a higher value of N r it even goes down to about 
10 -13 . As expected the numerical errors are most elevated 
at the boundaries of the different subdomains where the 
piecewise solutions of adjacent shells are matched. The 
GRV2 identity which had already proved its usefulness in 
the study of axisymmetric stationary neutron stars was 
applied to estimate the total error of the numerically com- 
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puted Kerr spacetime. The error estimator |1 — A| turned 
out to be closely related to the numerical errors derived 
above from comparison with the analytical solution and 
confirmed in this highly relativistic problem to be a sensi- 
tive indicator of the global numerical accuracy. 

6 Conclusion 

We have presented the application of spectral methods 
to several problems of numerical relativity. In each case 
they proved to be a highly valuable tool which lead to 
results typically several orders of magnitude more accu- 
rate than corresponding codes based on finite difference 
schemes. Especially in sphericallike coordinates the advan- 
tages of a spectral method which allows a rigorous treat- 
ment of the associated regularity conditions, while improv- 
ing the efficiency of the code at the same time, are remark- 
able. Particularly important properties for our problems 
are the negligible numerical viscosity in temporal evolu- 
tion schemes which enabled us to capture subtle details in 
the time-dependence of evolved variables as observed for 
equilibrium configurations of neutron stars in Sec. 3.3, as 



well as the very natural treatment of boundary conditions 
and the efficient solution of elliptic equations which is a 
frequently encountered task in our investigations. Our so 
far very positive experiences with spectral methods give 
us confidence to dispose of the appropriate numerical tool 
to tackle the exciting problem of black hole formation by 
3D-gravitational collapse of neutron stars. 
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